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Part I 



Introduction 



1. Physical background 



1.1. What is the most accurate number in physics? 

What is the most accurately known number in physics? It can be argued that theoretically 
this is the ratio between the perimeter and diameter of a circular disk, also known as vr, 
quite recently computed to about 10^^ decimals |,1|. This precision is not quite matched by 
experimental observations. It can also be criticized for not taking into account the discrete 
and quantum nature of matter, or — at this accuracy — even the tiny non-euclidean nature 
of surrounding space. 

Experimentally the best accuracy is probably what can be obtained by use of optical 
frequency combs [21, currently with a relative accuracy of a few parts in 10^"^ (i.e. 17 
decimals). This is expected to improve by a few orders of magnitude during the next 
decades, cf. Fig 2 of reference |2] Hence, there are situations where it makes sense to 
compute physical quantities to 20 decimals precision or better — provided that the physical 
model is known accurately enough. 

The latter is usually not the case. The physical quantity where experimental and theo- 
retical values are in best agreement is probably the magnetic moment of the electron. This 
quantity is measured to about 13 decimals l3j, with the computation of the theoretical 
contributions from Quantum Electrodynamics (QED) recently completed to fifth order in 
the fine structure constant a [4]. This leads to agreement between experiment and theory 
to about 13 decimals without adjustable parameters. Also the QED contributions to the 
muon magnetic moment has been computed to fifth order in a p], but for this quantity 
the theoretical contributions from other sources (like hadronic contributions to vacuum 
polarization) are larger, and the experimental uncertainty is also larger. 

1.2. The role of ordinary differential equations in physics 

Most of physics can be described locally in space and time, hence mathematically by 
differential equations. Differential equations therefore form a central part of theoretical 
physics, both on the elementary and advanced level. In most cases the relevant equations 
are partial. But by symmetry reductions, or more systematically for linear equations 
by the method of separation of variables, they can be reduced to ordinary differential 
equations. This greatly increases the prospects of finding solutions, and of understanding 
the properties of such solutions. 

There is a powerful and quite complete method of solving ordinary linear homogeneous 
equations, starting with works by Fuchs l6j and Frobenius l7j. In this method the solution is 
expanded in a convergent (generalized) power series, a Frobenius series, around ordinary or 
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regular singular points of the equation in the complex plane. This method can be combined 
with analytic continuation to extend the solution beyond the radius of convergence of each 
power series. 

1.2.1. Schrodinger equation in one dimension 

The Schrodinger equation is a partial differential equation that governs the time evolution 
of quantum mechanical wave-functions ^(q,t), where q denotes the position(s) of the 
particle(s). It gives a good description of the quantum motion of non-relativistic particles 
interacting instantaneously with each other. In the one-dimensional single-particle case it 
reads 



2m dq 



2+yiQ) 



^(g, t) = ih- 



dt 



1.1 



This equation can be reduced to an ordinary differential equation by separation of variables. 
Assuming a solution in product form, 

^>{q,t)=^{q)T{t). 

and substituting into equation (1.1), one finds T{t) = e~*^*/", and 

2m 5(?^ 



.2 + ^(«) 



^{q)=Ei;{q). 
We assume the potential to be a low-order even polynomial, 

V{q) 



N 
n=0 



^1.2) 



For numerical computations one must use dimensionless quantities. We first introduce a 
dimensionless length x = q/X such that the Schrodinger equation becomes 



fl2 N-i 

n=0 



ip{x) = eip{x), 



:i.3) 



with a2^+2 = h'^/2mVN, Vn = 2mT4A2"+V^^ and e = 2mE}?lf?. Here the scaling 
coefficient A has been chosen to give unit coefficients in front of the d"^ /dx"^ and the x^^ 
terms. Other choices may sometimes be more convenient, in particular it is natural to 
generalize to the case where 

dx"^ dx'^ ' 

with s^ usually a small number. It is common to think of it as fi? /2m^ but since this 
quantity is not dimensionless it is not a true small parameter of the equation. 



The solutions of equation (1.3) be expanded in a power series in x . The radius of 



convergence of this power series is infinite. One may make an analytic continuation to 
expand the solution around another point xq, but this will destroy the explicit parity 



symmetry of the problem. The latter leads to a doubling of expansion coefficients in 
the power series, and a significant increase in the number of coefficients describing the 
polynomial potential. The advantage for a numerical evaluation is that the power series 
may converge faster, and with a smaller loss of precision due to roundoff errors. 



1.2.2. Schrodinger equation in higher dimensions 

The real world is not one-dimensional, but can in many situations be treated as quite 
symmetric. F.i., the Schrodinger equation for a D-dimensional rotationally symmetric 
system. 



2m 



V^ + V{q) 



^{q,t)=ih-^{q,t), 



allows for a separation of variables. 



[lA) 



:i.5) 



Here q= \q\, and 3^^ '(g) is the generalization of the spherical harmonics to D dimensions. 
They are independent of the length of q (i.e., scale invariant), hence q • V3^*-^^(q) = 0. The 
functions 

V^'\q)^q'y^'\q) 

are homogeneous polynomials of order i in the components of the vector q, which are also 
solutions of the Laplace equation 



^1.6) 



V2/3^W (q) = {V\') 3;W (q) + q' V' D^^ (q) = 0. 
It follows from these relations that 

V^ [3;W (q)] = - {q-' VV) y^'^ {q) = -^{i + D - 2) q-' 37^ (q) 



and 



V2 3^W (g) i;{q) = 3^^ (q) [i>"{q) + {D - l)g- V'(9) - i{£ + D - 2) q-'i;{q) 



Hence, the separation of variables (1.5) leads to the radial equation 



!}L(^ -P-l ^ i{i + D-2 ^ 

2m \dq'^ Q dq q^ ' 



ij{q)=Eij{q). 



[1.7) 



This equation has a regular singular point at g = 0, but it still has a generalized series 
solution which can be found by the Frobenius method [7]. 

Separations of variables can be applied to many other coordinate systems. F.i., in 
three dimensions the Schrodinger equation in zero potential can be separated in ellipsoidal 
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coordinates (1^1,^2,^3), related to Cartesian coordinates by 



a2(a2 - 62) 



V 62(52 _ ^2) 

ab 



, 6 > a > 6 > 6 > 6 > 0, (li 



plus 10 degenerate forms of these coordinates fs]. The separated equations have five regular 
singular points, at iba, ±6, and 00 |9|. The Schrodinger equation remains separable if we 
add a potential of the form 

{il-iim-iDiil-il) • ^^^ 

The degenerate forms often lead to situations where two or more regular singular points 
merge to irregular singular points (confluent singularities). 

1.2.3. The Fokker-Planck equation 

There are of course many other fields of physics where equations similar to the Schrodinger 



equation occur. One such example is the Fokker-Planck equation 10 11 describing the 
time evolution of a probability distribution p{r,t) of a diffusing particle in a force field 

F{r) = -VU{r), 

|- p{r, r) = ^ V^p{r, r) - V • [F{r)p{r, r)] . (1.10) 

By writing p{r) = e~^^'^' '^{r,t) we obtain a "Schrodinger equation", in imaginary time 
T = it, for "^{r,T), 

- |:^(r-,r) = -^ V^^{r,T) + V{r) ^{r,T), (1.11) 

with a potential V{r) = ^ [F{r)-F{r) + V--F(r)]. For this reason we have denoted the 
equations we study in Paper I |35| and Paper II 136] as Schrodinger type equations instead 
of just Schrodinger equations. 



2. The projects of this thesis 

Although the projects of this thesis are inspired by all the points in section |1.1[ they are 
nearest to the first one. In the first project [35] the ground state energy of the anharmonic 
oscillator was computed to the accuracy of one million decimal digits. There are certainly 
no experimental results which require such accuracy, or physical systems which is modeled 
to that accuracy by the given Hamiltonian. That project was an attempt to explore the 
borders of applicability of our method, using standard computers of the time. 

The practical applications of very-high-precision calculations is that they provide es- 
sentially exact results, and hence may be useful for testing theoretical conjectures or the 
accuracy of approximation methods. One may also think of practical problems where 
numerical accuracy to several tens of decimals may be useful. 

2.1. Properties of the used algorithms 

Most projects of this thesis are based on an algorithm for solving Schrodinger type equa- 
tions to very high precision. Despite the numerical character of this topic, the focus of 
all investigations has been to explore analytic questions. How do the required resources, 
like computer memory and CPU cycles, scale with the wanted precision? With standard 
numerical methods, based on discretization of the Schrodinger differential operator, the 



error e = W ^ scales like a low power of the discretization step 16 , 

e ~ (Ax)", 

usually with 1 < n < 6. Hence the number A^ of computational steps would grow 
exponentially with the wanted precision, 

N ~ 10^/'^. 

In contrast, in the algorithm used here the number of computational steps grows asymp- 
totically linearly with the wanted precision D as -D — )• oo, 

N ^ iVo + const X D. 

This algorithm works for Schrodinger type equations with polynomial potentials in one 
dimension, and can be extended to a much larger class of ordinary differential equations. 
Note that in many situations the offset Nq can be quite large, making it computationally 
very expensive to obtain the first few decimals of accuracy. Hence, our algorithm may 
not be suitable for cases where only standard double-precision is wanted (but there are 
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examples where it is a competitive method). And we have not yet seriously explored the 
possible improvements which may be made by appropriate use of analytic continuations. 
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Figure 2.1.: The most time-consuming algebraic operation in our solution algorithm is the 
multiplication of two high-precision numbers. This figure shows a measurement 
of how this time varies with precision L> on a standard linux workstation. The 
theoretical prediction is that it should grow like a power D^'^ for low values 
of D. As can be seen this is followed quite accurately for precision values in 
the range 100 < D < 10^. Eventually the theoretical growth should reduce 
to a rate D \ogD log log D; this does not seem to be achievable in practice 
on our computers. The discontinuous behavior at low D occurs because the 
real computational precision increased in steps of 64 bits, corresponding to 19+ 
decimal digits. 



It should be kept in mind that each computational step will require more time as one 
increases precision. In our method the multiplication of two high-precision numbers is the 
most time-consuming operation. We have based all high-precision numerical computations 



on the CLN Class Library for Numbers 12 , compiled with the GNU Multiple Precision 



Arithmetic Library |13| . These hbraries use the Schonhage-Strassen algorithm 14 for 
multiplying numbers. With this algorithm the time T of multiplying two high precision 
numbers scale theoretically at a rate between D^'^ and D logD log log D with precision D. 
Hence the total time requirement of our high-precision method increases somewhat faster 
than D^ with the precision D. In practical computations one finds T ~ D^'^ for D < 2000, 



cf. figure 2.1 



Although it is of limited practical interest to compute physical quantities like eigen- 
energies to much more than 20 decimals precision, it may be the case that computation of 
the wave-functions in some regions of interest require very precise knowledge of the eigen- 
energy — otherwise it would be impossible to find solutions with the correct behavior. 
This is one practical motivation for developing methods to solve eigenvalue problems to 
otherwise ridiculously high precision. 

In paper II |36j it was demonstrated, and analyzed, that the very-high-precison wave- 
functions can be inserted into straightforward numerical integration routines to compute 
normalization integrals to comparable precision within acceptable time. 

We believe that these algorithms also can be used to compute functional determinants 
of one-dimensional differential operators to comparable precision. Since the functional 
determinant of a separated sum of differential operators is the product of the individ- 
ual determinants, the algorithms can also be used to compute functional determinants of 
higher-dimensional separable differential operators [15| . The algorithms can also be used 
to compute the resolvent of one-dimensional differential operators, but in this case an 
extension to higher-dimensional separable operators becomes more complicated since one 
must first compute the partition functions of the individual one-dimensional operators. 

2.2. Method for solving the differential equation 

Our method of solution is (perhaps disappointingly) simple: We just expand the solution 
ip{x) in a power series, 

i;{x) = x''Y. «-a;™, (2.1) 

m>0 

where o" = for the even parity solutions, and a = 1 for the odd parity solutions. For a 



given X the quantities A„i{x) = amX^"^^'^ can be generated recursively from equation (1.3), 

1 ^ 

A.,i(x) = (,^^2 + a)(2.n + l + a) E/"(-) ^"^-"(-)> (2-2) 



where 



uo — e, for n = 0, 
y^(x) = { Vnx'^'', forl<n<A^, (2.3) 

X , for n = N. 



This recursion is initialized with A_„(x) = Suqx'^ {n > 0). As is seen from equation (2.2) 
only the last A^ -|- 1 coefficients A^ix) need to be considered at any time while the sum 
is accumulated. This means that one only needs to store A^ + 1 coefficients Am and at 
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most A^ + 1 coefficients Vn, together with the accumulated sums for ip(x) and optionally 
ip'{x). Hence the total memory requirement is 2N + 4 high precision numbers, which is 
quite modest when the potential is a low order polynomial. 



In principle the sum in equation (2.f ) runs over infinitely many terms. But since equa- 



tion (1.3) has no singular points in the finite plane the sum will eventually converge very 



fast, and can be cut off at some finite value M. The value of J\f depends on the desired 



accuracy of the sum, the value of x, and of course the parameters in equation (2.1). 

The energy eigen- function V'(^) should become very small when \x\ becomes very large. 
Hence there will be large cancellations in the sum, and associated accumulation of round- 
off errors. Therefore, the interesting aspects of this method is not to code the recursion 



algorithm (2.2), but f.i. to estimate the proper choice of A^ and the numerical precision 
which must be employed in the computation. We have found that the WKB approximation 
combined with a Legendre transformation can be used for this analysis, see article IV (a) 
(381. 



2.3. Eigenvalue conditions 

The standard eigenvalue condition is that ■iIj{x) — )• as x — )• ±00. With parity symmetry 
(x — >• —x) it is sufficient to consider only one of these limits. I.e., we may start with 
solutions ^lJ„(x^,£) of desired parity, and use the single condition 

lim ipa(x,£) =0, (2.4) 

to determine the allowed eigenvalues e. 

The assumption that the potential V should be even is not necessarily physical, and not 
really necessary. In the general case one can expand the general solution in two linearly 
independent solutions ipa{x;£) and c/9f,(j;;e), 

ip{x; e) = Ca fa{x] e) + Cf, (pb{x; e), 

where Ca and C^ are coefficients to be determined. The eigenvalue condition is then that 
there should be a nontrivial solution of the equation 

lim ( '^'^(^'^^ "Pbix-^e) \( CcL , Q 



x^ooy ip^[-x;e) ipb{-x;e) J \ Cb 
The condition for this is that the system determinant must vanish, 

hmdetf ^';(^'"\ ^^(^'"V 
x^oo \^ ip^i^-x;e) ipb{-x;e) 

lim [ipa{x;e)(pb{-x;e) - ipai-x;e) ipb{x;e)] = 0. 



(2.5) 



This is in principle not different from the condition that a single function must vanish, 
but is in practice a nontrivial extension. The numerical implementation becomes more 
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complex, since we must solve an eigenvalue problem, and the program execution will be 
more demanding since one must compute and store four functions at each step instead of 
one. The assumption that the polynomial V should be of low order, i.e. that A'^ should be 
small, is mostly motivated by storage requirement. 

With our method of solution it is not possible to evaluate the wave- functions at, or very 



close to, X = cxD. We must therefore replace the condition (2.4) or (2.5) with boundary 
conditions at finite, but sufficiently large, xq. For each eigenvalue £ there exist a Robin 
boundary condition at xq, 



which is equivalent to (2.4). We do not know R{xo) exactly (it even depends weakly on 



the eigenvalue e), but it can be estimated from asymptotic analysis of equation (1.3) as 
xo — > oo. Define 

N-l 

Qix) = x^^+ Y, vnx^^'-e, 

n=0 

and let x be the largest solution of Q{x) = (i.e., the largest classical turning point). 
Since Q{x) becomes large when x ^ x we may choose 

R{xq) = \/Q{xq) + higher-order corrections (2-7) 

for a sufficiently large xq ^ x. This may often be further approximated by R{xq) = c«. 
I.e., a Dirichlet boundary condition, 'iP{xq) = 0, at xq. The most important property is 
that the eigenvalue is not very sensitive to the precise choice of R{xq). For x > x the 
solution can be approximated as 



il^ix) « K+ie) Q(x)~i/^ exp ( T dt^/oit)) + 

^Jx (2.8) 



K_ (e) Q(x)-i/4 exp (- J_^ dt^t)) , 



where the coefficients K± are expected to be of the same magnitude in general (both of 
order unity and slowly varying with e). 

The exact quantization condition, the Dirichlet boundary condition at xq, and the Robin 
boundary condition at xq become respectively, with Qq = Q{xq) and Q'q = Q'{xq), 

K+(e) = 0, (2.9a) 

K+{e) + K^{e) exp (-2 ^ dt^fQ{t)) = 0, (2.9b) 

K+{e) - \K-{e) Q^^'^Q'o exp (y-2 f' dt^t)) = 0. (2.9c) 



To compare solutions let Eexact be the solution of equation (2.9a), and expand 



K+(e) = -fir+ (Sexact) {e - eexact) + • • • ■ 
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Then the solutions of equations ( 2.9bp and (2.9c) become respectively 



S^Dirichlet — S-exact 
^Robin — £exact i 



(2.10a) 
(2.10b) 
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Figure 2.2.: The exact boundary condition, that |V'(a;)| — )■ as |2;| — ;• oo, cannot be realized in 
a numerical computation, but must be replaced by a Robin or Dirichlet boundary 
condition at some finite (large) value xq. This figure illustrates how xq must be 
chosen to maintain a described accuracy 10^^ of the computed eigenvalues, for 
the pure anharmonic oscillator, V{x) = x^, for the lowest eigenvalue and eigen- 
value number A^ = 50 000. Large values of xq is computationally challenging, but 
possible, even if one evaluates the wave-function directly by a series expansion 
around x = 0. 



We note that the Robin boundary condition improves the accuracy by a relative amount 
oQq Qq, which may correspond to a few decimal digits. This may be further improved by 
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systematically adding correction terms to the parameter R{xo), cf. equation (2.7). 
order of correction will provide a few extra decimals of accuracy. 



Each 



This may be a good 

approach if one wants the eigenvalue to standard precision (i.e. 15-16 decimal digits) only, 
but it makes little difference if one wants hundreds of decimals or more. Then one must 
make use of the fact that the factoijj 



exp 



Xo 



dt\jQ{t) 1 ~ exp 



1/2^ 



2x0 Qq 

N + l 



(2.11) 



vanishes exponentially fast with increasing xq. As an example consider the case that 



Q{x) = X — £. Then the quantization condition (2.9b) or (2.9c) leads to an error in the 
eigenvalue of order 

-ho^M^e^ (2.12) 



Ae 



~ e 



which clearly vanishes exponentially fast as xq increases. Figure [2^ displays how large one 
must select xq to obtain a desired accuracy of eigenvalue number A'^. As can be seen, this 
value is quite large for large A^, even at moderate accuracy. This may not be a serious 
obstacle if ip(x) is evaluated by analytic continuation, but it is quite time-consuming if one 
evaluates it by a direct power series expansion around the origin. 



^For this crude estimate of the integral we make a partial integration, using Q{x) = 0, and the approxi- 
2NQ{t). 



mation t^Q(t) 



Part II 



Mathematical background 
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3. Linear Ordinary Differential Equations 



In this part the mathematical backgromid used for the thesis projects is discussed. Some 
of the more technical parts is placed in the Appendices. 

A linear differential equation is any differential equation that can be written in the 
following form, 



P-^i^^ZTZ^. +Pn-l{z) 

or more symbolically, 



d 



n— 1 



+ 



-Po{z) 



m = 9iz), 



Cf = g. 



(3.1) 



(3.2) 



Equation (3.1) is said to be of order n, it is called linear because £ is a linear operator, 

C{cih + C2/2) = ci£/i + C2Cf2 (3.3) 

when ci and C2 do not depend on z, and ordinary because it only involves one independent 
variable as opposed to a partial differential equation. If g = the equation is said to be 
homogeneous, otherwise it is inhomogeneous. 

3.1. Expansions around ordinary and regular singular points 

We will assume the coefficients Pk{z) to be analytic (usually polynomials) in the region of 



interest. The homogeneous version of equation (3.1 ) is said to have a regular singular point 



at zq if some of the functions pk{z)/pn{z) has a (perhaps higher order) pole singularity at 
zo, but such that each 

ak = lim {z - zoT-''pk{z)/pn{z) (3.4) 



is finite. In generic cases the solutions can be expanded in Frobenius series around zq, 



f{z) = iz-zorY.fmi^-zor, 

m>0 



(3.5) 



provided cr is a solution of the indicial equation 



fc=0 



a 



(3.6) 



However, there are exceptional cases when equation (3.6) has multiple roots, or when some 



roots differ by integers, for which one must modify the series (3.5 ) with logarithmic factors 
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The general analysis becomes rather complicated, with many special cases to consider. For 
this reason the published code in Paper III |37| excludes all exceptional cases. 

For second order equations the exceptional cases can be reduced to only two possibilities. 
Since it was difficult to find the general recursion formulas for even these cases in the 
literature, in an explicit form suitable for coding, we have derived and published them in 
Paper IV(b) [39]. 



3.2. Second-order equations as first-order systems 

Many ordinary differential equations originate by separation of variables from partial dif- 
ferential equations involving the Laplace operator, and will therefore be of second order. 
They can often be transformed to the form 



d 



d 



K^) :7r2 + 9(^) :7r + '^(^) 



(iz2 



dz 



f{z) = 0, 



(3.7) 



where p{z), q{z), and r{z) are low-order polynomials in z, or the inhomogeneous version of 
such equations. This equation is well suited for expansion in a Frobenius series (cf. equa- 



tion (2.11) 



Hz) 



m>0 



Q»rn,Z 



since only a few coefficients Om-fc are required for computation of each next coefficient 
Om+i- Although the series is assured to converge up to the nearest singular point of 



equation (3.7), it may be convenient to evaluate f{z) indirectly through one or more points 



Zi by analytic continuation. Analytic continuation of functions which satisfy a second-order 
differential equation is rather simple to implement, since the function is fully specified by 
just two complex numbers, f{zi) and f'{zi), together with the differential equation. 



Hence we want to make coordinate transformations of equations like (3.7), and to im- 



plement robust algorithms for such transformations. This is simple as long as we limit 
ourselves to translations, since a translation of the independent variable, z = zo + u, only 



transforms equation (3.7) to 



au^ du 



f{u) = 0, 



where p(n) = p{zq + u), q{u) = q{zo + u), r{u) = r{zQ + u) are polynomials of the same 
order as the original ones. However, we would like to include the full group of Mobius 
transformations, 

au + 13 

These are the most general transformations of the Riemann sphere which do not introduce 



with aS — /37 = 1. 



(3.8) 



new singularities. It is not straightforward to describe a class of equations of the form (3.7) 
which are invariant under Mobius transformations, in particular to computers. For this 
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reason it seems preferable to reformulate (3.7) as a system of first-order equations. This 



should be done with some care. The perhaps most obvious choice 

yo{z) = f{z), yi{z) = f'{z), 
introduces an irregular singular point z = oo, since (with u = 1/z) 

-i-yo{u~^) = ^yi{u~^). 

au w^ 

This is unwanted unless z = oo already is an irregular singular point. Instead, if equa- 



tion (3.7) has singular points at z = zq and z = zi 7^ oo, the choice 
yo{z) = f{z), yi{z) = C{z- zo){z - zi)f'{z), 

will not introduce new singularities. One may choose the constant C 7^ freely. By setting 
C = —Zi and taking the limit zi — )• 00 one obtains 

yo{z) = f{z), yi{z) = {z- zo)f'{z). 

3.2.1. Example: Reformulation of the hypergeometric equation 

Consider the hypergeometric equation 

z{l - z)f"{z) + [c-{a + h + l)z\ f{z) - ah f{z) = 0. (3.9) 

This is known to have regular singular points at z = 0, 1, 00, and no other singular points. 
Hence, there are three possible pairs of singular points which may be used to define yi{z). 
We choose z = and 00, and a vector y with components 



yo = f, 2/1 = zf, 
and obtain the system of first-order equations 

d /yo\ / 1- z 



z{l-z) 



dz \yi 



ahz 1 — c+ {a + b)z 



(3.10) 



(3.11) 



The regular singular points of this equation, with the corresponding indices, can be ar- 
ranged according to the tableau 



1 00 

0a 

1 — c c — a — b — 1 b 



> . 



The first line of this pattern lists the positions of the regular singularities (here z = 0, 
z = 1, and z = 00), and the two entries below each position are the indices z^^ at that 
position. Here vi = 0, 1/2 = 1 — c at z = 0, ui = 0, 1^2 = c — a — b — 1 at z = 1, and vi = a, 



20 Linear Ordinary Differential Equations 

1/2 = b at z = oo. Note that the sum of all indices at all singular points is zero, not equal 
to one as in the 2nd order formulation. 

The general equation with regular singular points at 0, l,oo is obtained by considering 

f{z) =z''{l-zyy{z). Since 

z{l - z)-^zf'{l - zY = -z^(l - zY [/i(z - 1) + vz] , 
one finds 

The right-hand side is still a first-order matrix polynomial, in contrast to the second-order 
formulation where the polynomials increase in order when one generalizes the hypergeo- 
metric equation in the same manner. 



Finally, a Mobius transformation (3.8 1 which transforms the points 0, 1, cxd to mqi ^i, ^oo 
leads to the equation 

R\{u — Uo)(u — Ui){u — Uoo)^-/(w) 

au 
-vR'iiu-uo) + ^.{ui-u) Mi-U \ fiu) (3 13) 

with 

R, = (^—), R2 = ('^^^^), i?3 = (^^i^^). (3.14) 

V Uoo - no / Kuq- Uoo^ \Uo- Uoo^ 

The important feature here is that the factor multiplying ^ is a third-order polynomial in 
u, and that right-hand side involves a first-order matrix polynomial in u. This structure 
is stable under Mobius transformations, and transformations 



3.3. System of first-order equations 

Consider a system of first-order equations. 



p{z)^^y{z) = A{z)y{z), (3.15) 

where y{z) is a A'-component vector, p(z) is an ordinary polynomial of order N, and A(z) 
is a K X K matrix polynomial, 

N N-1 

p{z) = Y^p^z^ = C X{{z-zu), (3.16) 

fc=0 fc=0 

A{z) = Y.Akz\ (3.17) 

A:>0 
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3.3.1. Expansion around an ordinary point 

If Po 7^ the solution can be expanded in an ordinary power series, 

y{z) = J2 a-mz"". 



m>0 



We insert the series into equation ( |3.15 ) and introduce matrices 
In terms of these one finds the recursion formula 



(3.18) 



(3.19) 



(3.20) 



where ao can be chosen freely, and a-n = for n = 1, 2, • • • . According to general theory 



the series (3.18) will converge at least to the closest zero of p{z). I.e., the radius R of 

R>mm\zk\. (3.21) 

k 



convergence satisfies 



3.3.2. Expansion around a regular singular point 



If Po = but pi ^ the point z = is a regular singular point for equation (3.15). The 



solution can be found by use of the Frobenius method. We assume a solution of the form 



yiz;!^) = XI "« 



z 



m+u 



(3.22) 



m>0 



and find that equation (3.15) implies 



p{z] 



dz 



Ajy{z]v)= X I ^Mk{m + v)am-k 

m>0 \fc>0 



m-\-u 



(3.23) 



where the coefficients a__„ = for n = 1, 2, • • • . The coefficient of each power z'^^^ must 
vanish. For m = this implies 

M^{v)aQ = 0, (3.24) 



which has a nontrivial solution only when the indicial equation, 

detA4o(i^) = 0, 



(3.25) 



is fulfilled. This ii'* -order algebraic equation has K solutions Vs (counting multiplicities). 
There is at least one right eigenvector ao(z^s) for each distinct index Vg, and equally many 
left eigenvectors ao(z^s)- Each index Ug corresponds to an eigenvalue A^ = VsPi of the 
matrix Aq. 
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3.3.3. Distinct indices with non-integer differences 

Assume first that all indices Vs ai'e distinct, and that the difference between any two of 
them is non-integer. Then higher-order (vector) coefficients can be computed recursively 
as 

am+i{i') = -Mo{m + 1 + u)-^ ^ Mk{m + 1 + u) a„+i_fc(i/), (3.26) 

fc>i 

for m = 1, 2, • • • , where the coefficients a_„ = for n = 1, 2, • • • . The recursion is solvable 
at each step since all matrices A^o('7i + 1 + 1^) are invertible by assumption. (If one of them 



were not, the corresponding m + 1 + u would also be a solution of (3.25), contrary to the 
assumption.) 

3.3.4. Distinct indices; one pair with integer difference 

Assume next that all indices Vs are distinct, that the difference between two of them is 
integer, f2 = z^i -|- ^ with ^ > 0, and that all other possible differences are non-integer. We 
must then make a more general solution ansatz for the z^i-solution, 

y{z; ui) =Y.am z^+''' + Y.h„, z™+-i log(z). (3.27) 



m>0 mX. 



Now equation (3.15) implies 



d 



p{z)—-Ajy{z;vi)= ^ ^ 7Wfc(m + i/i) a„_fe z™+^i 

m>Ofc>0 (3 28) 

+ E E Pk-^i^m-k ^"+"^ + Mk{m + u^) b^_k ^'"^"^ log(z) = 0. 

m>ek>0 

The coefficients of each of the terms z'^^'^ and 2"^+^^ ^og{z) must vanish. For m = this 
implies 

Mo{ui)ao = 0, (3.29) 

as before. For m = i this implies 

Mo{i + iyi)be = 0, (3.30) 

which has a nontrivial solution since i + i^i = z^2 is also an index by assumption. Hence 
bi must be chosen as a right eigenvector corresponding to the index 1^2 ■ There is a corre- 
sponding left eigenvector bi. For m = £ we must further have 

Mo{i + i^i) ae = -pi be-Y, ^k{( + vi) be-k- (3.31) 

fc>i 

This is not solvable for a£ in general, since det A^o(^ + ^1) = 0. The solution criterion is 
that the right-hand side must be orthogonal to the left eigenvector b^, i.e. that 

Pih-be + Y^bi- Mk{i + 1^1) be-k = 0. (3.32) 

k>i 
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We can always choose the length of bi such that this equation is fulfilled. Then (3.31 ) can 



be solved for a^. The solution is not unique. We can add a right eigenvector b^ of arbitrary 
length to a£, and still have a solution. This corresponds to adding a solution proportional 
to 



y{z,U2) = X! f^rn{l^'i) 



m+U2 



m>0 



The remaining coefficients can then calculated by the recursion formulas 

am+ii^i) = -Mo{m + 1 + i/i)"^ Y^ Mk{m + 1 + i^i) a„+i_fe(i^i), (3.33) 



fc>i 



for m + 1 = 1, 2, • • • ,i — 1, and 



1^1 



(3.34) 



am + l(!^l) = -Moim+ 1 + !/l) ^ 2j-'^fc(™-+ 1 + J'l)am+l-fc(!^l) +Pft + 1 bm + l-fc(!^l), 

fc>l 
for m + 1 = £ + !,£ + 2,---. 

3.3.5. One double degenerate index 

Assume next that an index vi is doubly degenerate, that all the others are distinct, and that 
all index differences are non-integer. If there are two linearly independent right eigenvectors 



<im(i^i,?^) corresponding to the index i/i, then the solution ansatz (3.22) still works: 



m>0 r=l 



m+v\ 



(3.35) 



If there is only one eigenvector one again makes an ansatz with a logarithmic term 



y{z-vi) = Y am{vi)z^+-^+b^{u{)z^+-^ log(z). 



(3.36) 



m>0 



Equation (3.15) implies that 



MQim + vi) bm + Y -^kim + ui) b^-k = 0, 



k>l 



(3.37) 



Mo{m + ui) am + Pibm + ^ A^fc(m + i/i) a^.^ +pfc+i&m-fc = 0, 



fc>i 



for m = 0, 1, • • • , with 6_„ = a_„ = for n = 1, 2, • • • . For m = this implies that 

7Wo(z^i)&o = 0, (3.38) 

i.e. that &o rnust be a right eigenvector corresponding to the index z^i, and 

Mo{ui) a^, = -pihQ. (3.39) 

Even though detA^o('^i) = 0, this equation does have a solution as explained in subsec- 



tion 5.1 It is not a unique solution, because we may add a right eigenvector feo of arbitrary 



length to tto, and still have a solution. 
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3.4. Construction of the resolvent 

For a real independent variable x an integral expression for a solution of the inhomoge- 
nous equation can be written down if n independent solutions {yk{x) \k = 1, . . .n} oi the 
homogeneous equation are known. We first find a solution G{x, x') of 



P^^^'^TZ^ +^«-i(^)x:7^ + • • • +^o(^) 



dx"^ 



dx'^' 



G{x,x) = 5{x — x). 



(3.40) 



A solution of equation (3.1) can then be expressed as 

f{x) = / G{x,x')g{x')dx'. 



(3.41) 



Since G{x, x') is a solution of the homogeneous equation when x 7^ x' we must have 



G{x,x') = ^CkVkix), 
fc=i 



(3.42) 



where the coefficients Ck depend on x' , and is different for x > x' (denoted c^) and x < x' 
(denoted c^). Only the difference Ac^ = c^ — c^ contributes to the inhomogeneous 
solution. The function G(x, x') and its first n — 2 derivatives with respec t to x (denoted 
G^'^''(x,x')) must be continuous at x = x'. Further, by dividing equation (3.40) by Pn{x), 
and integrating it from x = x' — e to x = x' + e we obtain one more condition. Altogether 



lim [G(*^)(x' + e,x')-G('=^(x'-e,x')l =0, for A: = 0, • • • , n - 2, 

e— >0+ '- -• 



and 



(3.43) 



lim [G("-i)(x' + e,x')-G("-^)(x'-e,x')] = Pn{x')-\ 



This can be written in matrix form, 



yi 



Vi 



(1) 



y2 



(n-l) (n-1) 





(3.44) 



where all quantities are evaluated at the point x'. According to Cramer's rule the solution 
of this equation can be expressed in terms of determinants. For n = 2 one finds 



AcA ^ 1 (-y2{x')\ 

AC2J W{yuy2){x')p2{x') V 2/1(^0 J ' 



(3.45) 



where W{yi, y2){x') = \yi{x')y2{x')—y'i{x')y2{x')] is the Wronski determinant. By choosing 
c^ = C2 = Q this gives for n = 2, 



G(x,x') 



yi{x<)y2{x>) 
W{yi,y2){x')p2{x'y 



(3.46) 
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where x< = min(a;, x') and x> = niax(a;, x'). 

G{x,x') is the kernel of the integral operator C^^. This method is often referred to as 
variation of parameters . It can be extended from equations formulated on the real line to 
equations formulated on well-behaved curves in the complex plane. 

For equations formulated in regions of the complex plane one should instead search for 
a solution to the problem 

CG{z,z') = —^, (3.47) 

[z z) 

in terms of which 

r dz' ~ 
f{z) = J^—G{z,z')g{z'), (3.48) 

will solve Cf{z) = g{z) when C is a suitable curve in the z'-plane, encircling the point z 
once in the positive (anticlockwise) direction. 



4. The WKB approximation 



It is almost 100 years since Niels Bohr made his first formulations of a quantum the- 
ory of matter 17 18). Subsequent developments by him and others, including William 



Wilson [19l and Arnold Sommerfelt [20|, completed the "old quantum theory", and led to 
the formulation of the Bohr-Sommerfelt or Sommerfelt- Wilson or Bohr-Sommerfelt- Wilson 
quantization rules. With the advent of quantum mechanics and the Schrodinger equation 



these rules can be derived by the WKB approximation 21-241 to some extent. 

In this thesis we have used the WKB method for many purposes: (i) To make a priori 
estimates of the wave- functions before computing normalization integrals, as done in Paper 
II [36j , (ii) to estimate the magnitude of coefficients in the Frobenius series by use of the 
leading order WKB approximation, as done in Paper III [37] , Paper IV(a) |38| and Paper 
IV(b) |39| and (iii) to compare the very-high-precision numerical solutions against higher 
order WKB results, as done in Paper V |40|. 

The time-independent Schrodinger equation usually has the form 



"ip (x) = Q{x)ilj{x) 



(4.1) 



witlrj Q{x) = V{x) — E. We will mostly consider cases with the boundary condition 
lima;^-l-oo i^{x) = 0, and where V{x) is a polynomial in x. 

The standard WKB formulas can be found in most textbooks on Quantum Mechanics, 
A more thorough discussion is given by Bender and Orszag 12^. Then two 



f.i. 25 26 



leading order WKB solutions of equation (4.1) are 



i^±{x) = Q{x)~^'^ exp (±^ y^ ^Q(t) dt 
To satisfy the boundary conditions to this order the quantization condition 



(4.2) 



(4.3) 



0, 



must be fulfilled in the case of a situation with two turning points x± at which Q{x± 
so that Q{x) < for x- < x < x^. Here A^ = 0, 1, • • • . 

There are two extensions of these formulas which may be less known. These are (i) the 



Langer correction to (4.2) and (4.3) near a (regular) singular point, and (ii) a more general 
quantization condition derived by Dunham l28J which makes higher order WKB corrections 



to (4.3) quite straightforward to compute. These extensions are discussed in this chapter. 



^Note that this notation differs from the corresponding one in paper HI 
IV(b) [39] 
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paper IV(a) |38|, and paper 
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28 The WKB approximation 



4.1. The Langer correction 



Consider a generalization of equation (4.1), 



dx^ 



1 — Uj^ — V- d Vj^V- 
X dx x^ 



'iIj{x) = Q{x) 'iIj(x). 



(4.4) 



Such equations may f.i. arise as radial equations of rotation symmetric problems. The 
difference is that we now have a boundary co ndit ion at x = 0. By writing x = e" the point 
X = is transformed to u = — oo. Equation (4.4) becomes, with il^{x) = Q2'y^+~^^-i'^ '^{u), 



du 



2 ii^+-^-y 



m{u) =e^"Q(e")^(n). 



(4.5) 



This equation can now be solved by the WKB method, and transformed back to the x- 
variable. The results are that 



V;±(x)«x-±(Q(0)/Q(x))'/' exp(±^" 



^/Qit)-\[W) 



f) 



(4.6) 



Here Q{x) = 46 (i^+ — V-) + x Q{x). The normalization has been chosen so that 

il)±{x) ~ x"^ as X —7- 0^. 
The Langer corrected quantization condition becomes 



£ Jx- r 



N+\]t^. 



(4.7) 



where Q{x±) = and Q{x±) < for x_ < x < x+. This corresponds to the classically 
allowed region, with x± being the classical turning points. For the radial Schrodinger 
equation in 3 dimensions. 






(4.8) 



with e^ = f}^ /2m^ we have i/+ = £, i/_ = — (£ + 1). The quantization condition becomes 

- r* Je - V{r) -e^ii+ J)2 r-2 dr = (iV + J)7r. (4.9) 

In this case the Langer correction is a modification of the "centrifugal potential", 

v2 



'+1) 



+ 1 



With this correction the WKB quantization formulas for the hydrogen atom and the 3- 
dimensional rotation symmetric harmonic oscillator turns out to be exact. 
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4.2. Higher order WKB quantization condition 

4.2.1. Recursive calculation of higher order corrections 

Higher order corrections to the WKB approximated wave-function, 



ipix) = exp 



Ee"'5n(a 



n=0 



(4.10) 



can be found by substituting (4.10) into (4.1), and comparing terms order-by-order in e. 
We have the equations 



We find recursively. 



•^0 ~ Q^ 

n 

S'^_^ + ^SjSl^_j = 0, forn = l,2,.... 



5*0 — — vQ, 



n-l 



'"" 2VQ r ""' ^ § '^^'^"~^' 



(4.11a) 
(4.11b) 

(4.12a) 
(4.12b) 



The first terms of the recursion (4.12b) are 
lQ'{x) 



S[{x) 
S',{x) 



1 d 



4 Q{x) 2 dx 

5 Q'(a;)2 1 Q"{x 



logSo{x), 



32[Q(a;)]5/2 8[Q(a;)]3/2' 
IbQ'ixf 9 Q'{x)Q"ix) 



+ 



64 g(a;)4 32 Q(j;) 







(4.13a) 






(4.13b) 


1 Q'"(x) 
16 Q(x)2 


1 d 5^(x) 
2dxS'o{x)' 


(4.13c) 



There is another solution obtained by changing the sign of all even terms S2 



2m- 



The solutions can also be expressed in terms of multivariate polynomials. Define the 
infinite-dimensional vector 



Q = {Qo, Qir-- ,Qk,---) = {Q{x), Q'ix), ■■■ , Q^'^Hx), ■■■). 
Then the general term can be written as 



j„ 



[4Q( 



l(l-3n)/2 



TniQ), 



(4.14) 



(4.15) 



where TniQ) is a homogeneous n'th order polynomial in the components of Q, with integer 
coefficients, and also homogeneous of n'th order in derivatives. I.e., it consists of all 
monomials of the form 



QTQl 



-,no/^ni 



' Q^'' • • • , with all nfc > 0, ^ Uk 

k 



n, 



k 



n. 
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This means that number of terms in Tn is equal to the number of partitions of n. The first 
terms are 



To = -1, 

T2 = WQl-8QoQ2, 

Ts = -120 Ql + IUQ0Q1Q2 - 32QIQ3, 

7i = 2210 Qf - 3536 Q0QIQ2 + 608 QIqI + 896 QIQiQs 



128 QIQ4. 



(4.16a) 
(4.16b) 
(4.16c) 
(4.16d) 
(4.16e) 



Equation (4.13) indicate that each of the odd terms can be written as the derivative 



of expressions involving the even terms, and hence can be integrated explicitly. This is 
the case in general. In the classically allowed region, E > V{x) or Q{x) > 0, all even 
terms 82^ are imaginary and all odd terms S'2„_,_]^ are real. I.e, one can write the sum 



S' 



J2'n=0 ^"'^n ^S ^ 



Si, + iS[, with S'j, = E„^=o e'™+'5^ 



2m- 



_^ and Sj 






m=0' 



Jim 



CI 

Or, 



2m 



both real. They satisfy the equation 



{S'^^\S'jf^e{S'i-r\S'i) = -Q. 



From the imaginary part of this equation we find 



S'^ 



€ d 



log 5*1 



2dx 
which implies that (for rn, > 0) 



'2dx 



log 5^ 



' "^ log(l+£e2™^2m 



2dx 



S' 



(4.17) 



(4.18) 



On 



2m+l 



^ '^ log(l + f;e2'=^2fc 



2dx 



k=l 



nl 



(4.19) 



Order e^™ coefficient 



Hence 5*2^ ,]^ is the derivative of a single-valued function when m > 0. The next example 



beyond equation (4.13c) is 



Sk 




2 5^^ 



(4.20) 



This relation is straightforward to verify with a computer algebra program, but the explicit 
expressions in terms of Q are too lengthy to write down. 



One should be aware that the expansion (4.10) will not converge towards the exact 
result in general. Consider a case with a non-constant Q < everywhere, so that the 
WKB solution describes a wave moving to (say) the right. The higher order corrections 
will modify the shape of this right-moving wave, but never generate a left-moving wave. 
However, the exact solution for a quantum particle moving over a potential barrier will 
usually contain an exponentially small back-scattered wave. Hence one should expect 
exponentially small corrections to the WKB-series considered above. 
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4.2.2. The Dunham formula 



We now return to the two-turning point eigenvalue problem (4.1), with a potential V{x 



which is assumed to be analytic in x. The two-turning point quantization condition ( |4.3| ) 
has been generalized to arbitrary order in e by Dunham p^ (apparently as part of a Ph.D 
thesis at Harvard, after which no published research by the author seems to exist). 



1 

2k 



n=0 



■ Sj^{z)dz = Ntt. 



(4.21) 



The above integral is a complex contour integral which encircles a branch cut between the 
two classical turning points x± on the real axis. The WKB expansion breaks down near the 
turning points, but by extending the expansion into the complex plane the turning points 
can be avoided. The quantization condition is obtained by requiring the wave-function to 



be single valued. The integral in (4.21) is finite because the contour encircles the turning 



points instead of passing through them. The quantization condition (4.3) is recovered by 



considering the first two terms of the expansion. The contribution from 5*0 becomes 



1 

2ie 



1 
2ie 



SQ{z)dz - 
and the contribution from 5i becomes 



\/Qiz)d 



I rx2 , 

z = - J-Q(x)dx, 

e Jxi 



^/^iw"-'"-^ 



^logQMd. = -i4ri 



vr 
'2' 



since the logarithmic integral encircles two simple zeros. Hence the Dunham quantization 
condition becomes 



^'1 m=l 



{z)dz 



N+-]TT. 



(4.22) 



4.2.3. Exactly solved cases 



The quantization condition (4.22) seems to depend only on Q{z) in the region near the 
branch cut from xi to X2, and therefore obviously cannot always be correct. We could 
modify the potential in a far-away region, thereby changing the exact eigenvalues, without 
changing the value of Q{z) in the region of integration. However, this is only possible with 



a non-analytic potential. The interesting question is whether (4.22) is exact or not for 



analytic potentials. To our knowledge this had proven to be true for all cases where the 



expansion in (4.22) can be carried out to all orders, and a comparison with exactly known 



solutions can be made 28-31 



This is known to be the case for 



1. the harmonic oscillator. 



Q{z) 



E, 
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where the integration contour of all correction terms in (4.22) can be deformed to 



infinity, leading to the conclusion that all correction terms vanish. Hence, the con- 



dition (4.22) reduces to the standard first order WKB results, which is known to 



reproduce the correct result. 
2. the Morse potential p2], 



Q{z) =e-'^ -e-' -E, 
where the same procedure can be carried out after a change integration variable, 



u = e-^ in (4.22). 



3. the radial equation of the hydrogen atom, 

Q{z) = az^'^ + bz''^ - E, 
where the same procedure can also be carried out after a change of integration vari- 



able, u = z~^, in (4.22). 



4. the radial equation of the rotation symmetric harmonic oscillator. 



az^ + hz~'^ -E. 



One may show that only a subset of the terms in the WKB expansion have a non-zero 
integral. These terms can be computed explicitly. 

5. the Poschl- Teller potential j33l. 



Q{z) = asech z — E. 

By introducing the parameter u = sinhz one may show that only a subset of the 
terms in the WKB expansion have a non-zero integral. These terms can be computed 



explicity 29 



All cases have the common property that they have only a single branch cut in the full 
complex plane of the final integration variable. 

4.2.4. Polynomial potentials 



Now restrict to the case that Q{x) is a polynomial of order K. We observe from (4.13a- 
4.13c) that S^ has the form 



Snix) 



1 



[g(3.)](3n-l)/2 



Pn(x), 



(4.23) 



where P„ is a polynomial of order {K — l)n. By substituting this ansatz into (4.12b) we 
verify that it is correct, and find the recursion relation 



1 11" 

P„+i = -(3n-l)Q'P„--QP; + -5]P,P„+i_„ 



(4.24) 



i=i 



with Po = —1- This gives 
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Pi = -^Q\ 

32 8 



64^ 32 



QQ'Q" + ^Q^ Q'", 



(4.25) 

(4.26) 

(4.27) 

is the derivative of a single-valued 



1 



and so on. Since every odd S2mj^i,{m = 1,2,- 

function it does not contribute to the quantization condition (4.21 ). Thus, (4.21 ) simplifies 

to a sum over even-numbered terms only, 



1 
2k 



/OO -j 

^62-5;^(z)ciz=(iV+-)vr. 

m,=0 



Further, we may subtract any total derivative of the form 



d R2m{x) 

dx [Q(x)]"'" 



1 



-^ [Q{x)R2^{x) - am Q'{x) R2m{x)\ 



(4.28) 



(4.29) 



(with R2m a single-valued function) from S2m without changing the value of the contour 
integral. By choosing a^ = 3m — |, and i?2m a polynomial of order {K — l)(2m — 1), this 
can be used to replace P2mi of order {K — 1)2?tt,, by a polynomial P2m of order {K — 2) or 
less. 



Part III 



Appendices 
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5. Differential Equations 



In this appendix we provide a fairly complete analysis of the general expansion of solutions 
to a general first order homogeneous linear matrix differential equations at a regular sin- 
gular point, for the cases when higher order logarithmic terms occur. The derivations are 
done for the purpose of later numerical implementations, but we have not yet done such 
implementations. 

5.1. Matrices with fewer eigenvectors than eigenvalues 

An arbitrary D x D matrix fC cannot be completely diagonalized in general. Although the 
polynomial eigenvalue equation, 

det (/C - /x) = 0, (5.1) 

always has D solutions {^u^ \d = 1, . . . , D} counting multiplicities, a d-iold degenerate dis- 
tinct root /i may have dfj, < d linearly independent eigenvectors (with d^ > 1). I.e., the 
geometric multiplicity d^ of an eigenvalue /i may be lower than its algebraic multiplicity d. 
However, any D x D matrix /C can be brought to Jordan normal form |34|. I.e., if /C has 
d linearly independent eigenvectors it can be similarity transformed to a block diagonal 
form, 

Ji 

J = S-^1CS={ \ ■.. : 1 , (5.2) 

••• Jd 

where each block J„ is a d„ x dn bidiagonal matrix of form 



Jn=\ " ^"^ ■ M. (5.3) 





There may be more than one block for each distinct eigenvalue. Now observe that there is 
a sequence of column vectors, 

e(^) = (0,--- , ^ ,••• ,0)^, k = l,--- ,dn, 

position k 
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(where stands for transposition) such that 



(J„-/i„)eW =6(^^-1), for A; = 2 
( J„ - ^ln) e(i) = 0. 



) ) "'n 



There is also a corresponding sequence of row vectors, 

e^'^) = (0, • • • , 1 ,---,0), k = l,...,dn, 



such that 



position dn + l—k 

e(^) {Jn - /in) = e('=-i) , for k = 2,- ■■ ,dn 

e(l) {Jn - fin) = 0. 



(5.4) 
(5.5) 



(5.6) 

(5.7) 



Equation (5.5) means that e^"*^^ is a right eigenvector of J„, while equation (5.7) means 



that e^^-* is a left eigenvector of Jn, both with eigenvalue /i„. There is a (perhaps unusual) 
orthonormality relation for these vectors. Define the "backward identity matrix" E as 



E 



(5.8) 



Then we have the relation 



eW£eW=A 



ke- 



(5.9) 



One may extend the column vectors e'^ to D-dimensional column vectors g:('^'"'), and 
the row vectors e^''' to D-dimensional row vectors eC^'") by inserting them into the n'th 
block, with zeros in all other blocks. They satisfy the relations 



( J - /in) e^'^'"^ = e^''-''"^ for A; = 2,, 
(J-/z„)e(i'")=0, 



,dr, 



(5.10) 



(5.11) 



for n = 1, 2, • • • d. These relations, combined with equation (5.2), means that we can find 
a solution to the equations 

(/C - fin) a;(i) = 0, 



(/C - //„) a;(2) = C2 a;W 



(5.12) 



(/C - fin) a;^'^") = Cd„ X 



,(dn-l) 
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in terms of the vectors Se^^'"'\ We make the ansatz x^^'' = NkSe^^''^^ + PkSe^^'"'\ and 
find 

™(1) _ AT Cp(l-") 

k 

iE(fc) = iVi n Cj 5e(^'"), for /fc = 2, . . . , d„, 

i=2 

where A^i is a free parameter. 

5.2. General case 

In the general case the initiahzation step constitutes of finding a complete set of solutions 
to the equation 

(piz^ - Ao) a(z) = 0. (5.14) 

5.2.1. The Jordan normal form 

Aq can be brought to Jordan normal form by a similarity transform |34] , 

Ao = SJS-\ (5.15) 

with J' block diagonal 

/J. ••• OX 

\0 ••• J J 
Here each block J„ is a d„ x dn bidiagonal matrix of form 

'A„ 1 ••• 

Jn=\ ° ^" / " 1 • (5.17) 

••• 

This representation is unique up to permutation of the blocks. There may be more than 
one block for each distinct eigenvalue. Introduce Kn = Jn — ^n, and observe that (we let 
-^ denote transposition) 







(5.18) 
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This means that the matrix 



-T 



KnK^ =diag(l,l,...,l,0) 



(5.19) 



is a projection onto the subspace orthogonal to the left eigenvector e„ = (0, . . . , 0, 1) of 
Kn- For a coordinate invariant description, one notes by taking the scalar product of the 
equation 

Knf = g (5.20) 



with the left eigenvector e„, that 



BnKnf = = Bn-g, 



since e^Kn = 0. Hence, equation (5.20) has a solution only if the solubility condition 

en-g = (5.21) 

is fulfilled. Further, when there is a solution / it cannot be unique; we may always add a 



term aSn to /, since Kn Sn = 0. When equation (5.21) holds the general solution is 

f = Klg + aen, (5.22) 

where the coefficient a can be chosen freely. 

5.2.2. Solution of a homogeneous block equation 

The initialization step consists of finding the general solution of the homogeneous equation 

d 



PlZ- Jn) (p{z) = 0. 



(5.23) 



We make the ansatz 



^{z)= J2 ¥'fe^Mog(z)'^"-l-^ (5.24) 

fc=0 
where piv = An is the eigenvalue of Aq corresponding to the n'th block. Inserted into 



equation (5.23) this leads to the conditions 

Knifk = Pi (dn-k) ipk^i, for /c = 0, . . . ,(i„-l, 

where (^_i = 0. This means that (po must be an eigenvector of Kn, 

ipo = aoen = ao e^^^ 

The next dn — 1 equations can be solved iteratively, 

k 

(Pk = pi{dn-k)Klipk~i + ak e„ = Yl «J Pj ^n~^\ 

3=0 



(5.25) 
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for k = 1,- ■ ■ ,dn — l- Here 



eW = (Kj)e„, and P^ = l[ p^idn-i) with P|^ = 1. 

e=j+i 



(5.26) 



Hence we have found a (i^-diniensional space of solutions to equation (5.23) 

dn—l dn—l 



d„-l 



j=0 k=j 






z)"" ' "^ = ^ aj¥J^-"(z) 



(5.27) 



when pii' is an eigenvalue of J„. with algebraic multiplicity dn and geometric multiplicity 
1. 

5.2.3. Solution of an inhomogeneous block equation (regular case) 

The recursion step essentially consists of solving an inhomogeneous equation like 



PlZ— - Jn) (f{z) = J2 Xk Z^" \0g{zY ^•. 

"-^ fc=0 



(5.28) 



Assume that pi/i is not an eigenvalue of Jn-, so that pi^ — Jn is invertible. We make the 

solution ansatz 

e 

\i-k 



viz) = J2'Pkz''\og{ 
k=0 



Inserted into equation (5.28) this leads to the conditions 



(pi/x - Jn)<Pk = Xk-Piii + 'i--k) (pk^i with ¥?_! = 0. 
These can be solved recursively as 

Vk = {PlIJ' - JnV^ [Xk -Pl{(. + l-k) ifk-l] ■ 



(5.29) 



(5.30) 



(5.31) 



5.2.4. Solution of an inhomogeneous block equation (singular case) 

Assume that piu = A„ is the unique eigenvalue of J„, i.e. {piu — Jn) = —Kn- To solve 



equation (5.28) we make the ansatz 



dn+e. 



fi^) = H 'Pkz" \og{ 



\dn+e~k 



k=0 



which inserted leads to the condition 

Kn^Pk+i = -Xk+i +Pi idn+i-k)(pk. 



(5.32) 



(5.33) 



42 Differential Equations 



Here one should interpret ^-\ = 0, and Xfc = for A: = 0, 1 • • • , dn — 1- Since BnKn = 0, 



equation (5.33) can only have a solution if the right-hand side is orthogonal to e„. I.e., 



en • Xk+l = Pi (dn+i-k) Bn ■ (fk- 



(5.34) 



This determines the last component of all but one of the (/3fc's, 



^k 



c<*--»4. 



-' + ^1". 



for k = 0,l,--- ,dn+i-l. 



(5.35) 



Here 



C 



(dn-l) 



(e„ • Xk+i) 



(en ^ • Xk+l) 



pi{dn+£-k) pi{d„+i-k)' 



(5.36) 



and (pf^ denotes the restriction of cpk to the space orthogonal to e„. I.e. the first dn — 1 



,-.(1) 



components of ipk- We define ^^ in the same manner. Insert the partial solution (5.35) 

frfn-2) 



into equation (5.33), and use that K^en 



(rfn-l) 



,:;(!) 



. Since K^ifl. is also orthogonal to 
the space spanned by e^T" ^^ we can use this to determine the last component of all but 
two of the (pj^ s, 



Adn-2) 



n 



(1) 



c. 



(dn-2) (d„_2) 



+ if^^\ fork = 0,l,---,dn+e-2. 



(5.37) 



Here 



C, 



idn-2) 



Piidn+i-k) 



= (1) . ^(1) ' 
-n A.fc+1, 



+ c, 



(dn-l) 
fc+1 



(5.38) 



and (p^. denotes the restriction of (pk to the space orthogonal to e^ and Bn . I.e., the 

-(2) 

first dn — 2 components of cpk- We define xj^ in the same manner. Note that we cannot 



.(rfn-2) 



(dn-l) 



determine C^ "+e~i y^t, because we do not know C^ \t. ■ The solution process thus far is 



illustrated by the first two frames of figure 5.1 , for the case that 1 = 2 and d. 
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Figure 5.1.: This figure give a schematic description of how one can solve the recursion 
equations at a step m, when m + i^i is a multiple root of the indicial equation, 
and where one already have logarithmic terms in the expansion. 

~(r—l) 

This process can be continued. At step r the vector Kn^ff. is orthogonal to 

p(0) p(i) 

tin ) t^n 1 



-(r) 

, Cn , leading to 



-(r-l) 



c, 



(dn-^) Jdn'-r-) 



with 



C, 



id„-^) 



1 



pi{dn+i-k) 



+ (pt\ for A; = 0,l,...,d„+£-r, 






■'k+1 



(5.39) 



(5.40) 



After dn steps all consequences of the solubility condition (5.35) have been deduced, as 
indicated by the third frame of figure 5.1 The first (. + 1 vectors (fQ, . . . ,(p£ are com- 



pletely determined at this stage. To determine the last dn vectors completely we solve 



equation (5.33) in the forward direction, 

<fik = -Klxk +pi{dn + (- + l- k)K^ ipk^i + aken, (5.41) 

for k = i + 1, . . . ,i + dn. Each solution involves an arbitrary constant a^. This process is 



indicated by the last three frames in figure 5.1 



6. WKB Quantization of the Quartic 
Potential 



The contents of this appendix started out as a small example of the WKB quantization 
method to high orders. It soon grew in magnitude and resulted in Paper V 40 — and 



delayed the completion of my thesis with several weeks. It provides more details of the 
computations reported in Paper V, for the high order WKB analysis of the problem 



2 " I 4 



V'(x) = E7p(x). 



(6.1) 



One can set the parameter e = 1 without loss of generality, but it is useful for initial 
organization of the WKB expansion. 



6.1. High order WKB expansion 

As an example consider the case V{x) = x^. The first few polynomials P2m and P2m 
become 



n 


Pnix) 


Pn{x) 


2 
4 
6 


x^ + lEx^ 

14^12 ^ 3|3^^8 ^ 321^2 ^4 ^ |^3 

671:cl« + 3 223£;a;14 + 104 595^2^10 ^ 63 075^3^6 ^ 279^4^2 

4 S lb 2 


lEx-" 

77 p3 
1768 
61061 p4^2 
62 928-^ -^ 



The general pattern is that 



Pu{x) = {-lYE^'pf\ 



Pu+2{x) 



-lYE^'+^x^pf\ 



(6.2) 
(6.3) 



where the coefficients pf and pf are positive rational numbers. F.i., as can be seen from 

the table above, p^ = ^, p^ = j^, Pi = H^ ■ Thus we have to do two types of 
integrals 



r(e) 



r(o) 



2i 



2i 



1 



fc+l/2 ' 



dz, for k 



{z^ - E) 



— ^ dz, for k = U + 2. 
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We find, by writing z = E^''^u, and deforming the integral along the real axis, 



'-1 



'-EY^'dz 



E^l^ f duVl 



1 

2i 

2 M'2^ 3 M'2^ 



(6.4) 



where B(a, b) = T{a)T{b)/r{a + b) is the Beta function. The remaining integrals of interest 
cannot be deformed to convergent integrals along the real axis, but they are S-derivatives 
of integrals which can be deformed. For A; = we find, by writing z = E^'^^u, and deforming 
the integral along the real axis, 

f-i 



Ir 



(e) 



1 

2i 



1 



(.4 _ E) 



1/2 



dz = 2E-^I^ 



1 



^ 



u* 



<i«=^B(l,l)i.-/^ 



r(o) 



2i / f^4 _ K^l/2 L 



(z4 - E) 



vr 



du=iB(?,-)i?l/^ 
2 ^4 2 



By differentiating these relations A; times with respect to E we find 



Ke)^(_l)fcl(i)(i + l)---(i + ^-l) 



2(^)(^ + l) 



;Ufc-i) 



B(l l)i5;-i/4-fc 



(-l)\B(i i-fc)i?-i/4-fc, 



r(o) 



2-v4'2 

(-i)'i ^~~^!'^}^l}"'^~~'t^~.^^ B(|, i) i5;V4-fc 



(i)(i + l)---(i+fc-l) 



.2A2 

; 1 

^4' 2 



{-lf\^i.l\-k)E'l'-\ 



We introduce the quantity 



VT' 



■4^(i'i) 



1 r(i; 



-' r(f )^ 



7.764068 784... 



(6.5) 



(6.6) 



(6.7) 



in terms of which B(i, |) = -k p^l'^ and B(|, ^) = 4/)-i/4, Inserting ( |a2|a3| and ( |a5|a6 ) 
into (4.21) gives the quantization condition 



^pl/4 E^l^ - -p-y^ E'^l^ + ^ gf ^ ,4£-l^-(12£-3)/4 
oo 

+ E '?i°^e''^' i?-(i2^+3)/4 = (AT + 1 ) vr. 



^=1 



Here 



(e) 






^ ^ 2^^ (i)(i + l)...(l+6£-2) 
(-l)^2n^°-' ^~4A~4 T ^; • • • y-j -r u<, -r x; ^_^i^ 



)a + i)...a + Qi + i) 



(6.8) 

(6.9) 
(6.10) 



^2^V2 
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Borel transformed expansion coefficients B2irm) (scaled) 



"T 1 1 T 
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0.16 



0.14 



0.12 



0.10 



0.08 



0.06 



(m + 1)^/2 
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• • • Coefficients for odd series 
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Figure 6.1.: Rescaled form of the WKB expansion coefFicients rm of equation (6.12) for m 



0, . . . ,852. The index i/ = | is chosen as the simplest rational number close to 
the best fit, after which we find Or ~ 0.202 641 423 4 as the best fit to a sequence 
approaching a constant absolute value for large m. 



By introducing 



e = K e 



1,2 j^-3/2 



(6.11) 



with K = q p^'^ ~ 0.309 600 873 . . ., we can rewrite equation (6.8) as 



2£ „(e) 



with new coefficients r2£ = k q\ and r2i+i = k 



1 

,2^+1 „(o) 



(6.12) 



. The first few terms are 



1 



11 



n 



r2 



127r' 41472 

The further coefficients have the form 



P, r3 



4 697 p 



7464 960 vr' 



r4 



390 065 

8 026 324 992 



P 



r2E = (-1)^"^^ r2e /, r2e+i = (-1)^+^ r2e+i p'^/ir, 
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where fm are positive rational numbers. The coefficients Vm grow Uke m}?'a^ /{m + 1)^ in 
magnitude for large m. We have computed these coefficients up to m = 852 (corresponding 
to the 1 704* order of the WKB expansion) . Empirically they fit the cited behaviour quite 
well, with ar ^ 0.202 641423 4 and i^ = |, see figure 



6.1 



Obviously the sum r(e) = J2m=o ''»" ^"^ ™ equation (6.12 ) has zero radius of convergence. 
However, if one uses the integral formula 

/•oo rco 

^,2^^2(m+i)/ clxx'^e-"'' dyy"'e-'^y 

Jo Jo 

POO 

= 2a2("^+i)/ dC^"'Ko{2a^/^), with a = e'"^ (-f < </<< f ), (6.13) 

"'0 

and interchange summation and integration, one obtains an integral expression (Borel 
resummation) , 

POO roc °° 

r(e)=/ dxe""^ / dy e""*' V f ^ (xya^e)" , (6.14) 

Jo Jo r^n 

with 

rr 

Now the sum 



a 



m=0 

2(m+l) 



,2 m '^"'- (^-^^^ 

m\ oJIr 



{z)^Y.f^-.z^ (6.16) 

m=0 



converges for \z\ < 1. For a = 1 the function f[z) has singularities where z^ = —1, with 
the singular parts behaving like (1 + z^)'^' ^ near the singularities. In terms of the variable 
z^/(l + z^) this singularity is mapped to c«, and the full integration range is mapped to the 
interval [0, 1]. However, when one tries this substitution, in the hope that the (rewritten) 
sums for f{z) will converge over the full integration range, one discovers that there are 
additional singularites where z'^ ~ 4. Hence, to avoid integrating through a singularity, 
one must introduce the phase a (or equivalently integrate along a different direction in the 
complex plane). A convenient choice is a = e"^'^, or its complex conjugate. Actually, to 
assure a real result after analytic continuation of f{z) beyond the radius of convergence of 



the sum (6.16), one must take the average of these two choices. This amounts to taking 



the real part of the integral (6.14). 



After this choice we separate f(z) into four (infinite) sums. 



p=o e=o 



The function defined by each infinite sum is singular at z"^ = —1, z'^ ~ —16, and probably 
at infinitely many more points on the negative real z^-axis. Now rewrite 

J2 ^4^+P ^^^ = J2 ^4f+p (i^) , (6.17) 

e>o i>o 
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and use the computed coefficients r4e+p to find equally many coefficients r/j^i^p. By com- 
puting the sequence of coefficient ratios 



Pe 



(p) _ ?'4{^-l)+p 



r^e+p 



(6.18) 



we find empirically (by the ratio test) that the right hand sum of equation (6.17) converges 
for 1-2^/(1 + z'^)\ < 1, see figure 6.2 This completes the construction of a function 



r(e) = / dxe 
Jo 

roo 



OO /"OO 

ax 



dye °'^ f{xyar£) 



POO 

2 dCKo{2a^)f{Care) 
Jo 



(6.19) 



which has the same asymptotic expansion as the sum in (6.12), and where the integrand 



notably f{z), is computable to high precision over the full integration range. 

Coefficient ratios for the sums defining r(z) 



0.06 



0.05 



0.04 



0.03 



0.02 



0.01 



o.ao 



I 



p = 

p = 1 

p = 2 
p = 3 



M _ T _ ^4(^-l)+P _ 
P(^ J- — .". J- 




000 0.005 0.010 0.015 0.020 0.025 0.030 0.035 0.040 



Figure 6.2.: The convergence radius of the right hand sum in equation (6.17) is, according 
to the ratio test, given by the inverse ratio pj, of successive terms in this sum 
as ^ — )• OO. The computed ratios give convincing evidence that the convergence 
radius is unity. 
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6.2. Series solution of the quantization condition 



The quantization condition (|6.8|) can now be solved numerically by first initializing e 



next iterating the recursion 



A^ + 



r{ef 



(6.20) 



until convergence, and finally applying the relation (6.11) 



One may also proceed analytically by expressing e as a series in the (small) quantity 



5^{N+\r\ 



e = 5+Y,Sm5'' 



(6.21) 



m=2 



The coefficients Sm can be computed recursively. The first terms are 

1 



S2 = 2ri 



67r' 



S3 = 2r2 + br^ 



11 5 

P + 



20 736 



144 7r2' 



(6.22) 



Computation of the exact Sm, which are polynomials in p and vr^^ with rational coefficients, 
becomes too memory- and time-consuming beyond the first few tens. We have computed 
the sequence exactly up to ssg, and higher Sm with about 3 800 decimals accuracy. 



When the sequence of Sm is known one may use (6.11 ) to express i? as a series in 5, 

-2/3 



E^Ej, = K-^/^e'"U-^'^{l+Y. 



.i5" 



m>l 



^-2/3^4/3^-2/3(^1^^^^^™^ 



m>l 



^^-2/3^4/3 5-2/3 ^(^)_ 



(6.23) 
(6.24) 



The first few terms are 



h 

t2 

u 



1 

9^' 



11 



648 7r2 31104^' 
11 341 p 



8 7487r3 466 560 vr ' 

1 309 9 163 

+ 



P 



+ 



1 748 093 



5 038 848 7r4 25194 240 7r2 27088 846 848 



-,P ■ 



(6.25) 



In this way an expansion of the WKB-solution (4.10) to order 2N in the quantity e can be 



used to compute the expansion (6.23) to order N in the quantity {N + g) 



l\-2 
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Borel transformed expansion coefficients B2(tm) (scaled) 
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Figure 6.3.: Rescaled form of the expansion coefFicients tm in equation (6.23) for m 



0, ...,852. The index i/ = | is chosen as the simplest rational number close 
to the best fit, after which we find at ss 0.202 641423 4 as the best fit to a 
sequence approaching a constant absolute value for large m. Note the similarity 
with the r^-sequence. 



6.3. Extended Borel summation of the asymptotic series 



The sequence of t^ looks very similar to the sequence of r^, cf. figure 6.3 Hence one 
may use the same method to construct a convergent expression which reproduces the series 
expansion of t{6). 
We define 

Q,2(m+1) 



771! a"' 



and coefficients tm such that 



i{z)^ Y^tmZ'^ = Y,zPY.hl^ 



-P U+^4 



(6.27) 



m>0 



p=o e>o 
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One can use the previously computed coefficients tm to compute equall y man y coefficients 
im- With the chosen value of a = e"^'^ the sums over i in equation (6.27) converge for 



< z < oo, see figure 6.4 



Then, the expansion of the integral expression 



t{6) = Re\l dxe 



oo fOO 

ax 



dye "'y i{xyat5) 



(6.28) 



as a series in 6 reproduces the sum in equation (6.23). However, this particular integral 



may not be the best way to use the series expansion. 

Magnitude of coefficients tm (scaled) for the sums defining i{z) 
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Figure 6.4.: Scaled versions of the expansion coefficients in equation (6.27). They provide 
convincing evidence for convergence of the sums over the full integration range, 
< z < oo. 



A better approach is to regenerate the series expansion with an integral expression for 
the remainder. To this end write 

e-°^ = -a*-^e-°^ 
dx 
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in equation (6.28), and perform a partial integration, 



oo roo 

ax 



oo roo 

ax I J -ay 7(1) I 



dxe""^ / dye-°H{xyat6) = to + a*at6 / dxe""^ / dyye'^'y t'^'^xyatS). 
JO Jo Jo 



By repeating this process M times, and taking the real part, one finds 



M-l 
m=0 



(6.29) 



with 



Here 



{/■oo /"OO 



t^''\z) 



d^ 

d^ 



i{z) = ^ (m + M)(m + M - 1) • • • (m + l)tm+M z"' 



m>0 



^Ei 



(M) ^m 



m>0 



(6.30) 



(6.31) 



tion (6.27), 



From Ai known coefficients im one finds Ai — M coefficients im . Again rewrite, cf. equa- 

(6.32) 



tWfz) 



m>0 



■(M) ^m 



n-'Uf'.rUj 



p=0 f>0 



and use the known coefficients tm to compute equally many coefficients im ■ 



We finally insert the expansion (6.32) into (6.30) and perform the integral numerically. 
This gives tcorl{5) to a relative accuracy of about 10^"'^'^. As a consistency check we verify 
that t{5) is independent of M, at least for M- values around the point where tcorr(<5) 
is minimum. As can be seen qualitatively from figure 6.5 this works well for the lowest 
eigenvalues. But it also shows that the WKB-series does not reproduce the exact eigenval- 



ues, even when the correction term (6.30) is included. The quantitative results are shown 



numerically for the two lowest eigenvalues in tables 6.1 6.2 

The WKB-series shows a quite stable result when the correction term from Borel re- 
summation is added, with an uncertainty much smaller than the distance to the exact 



result. This can be seen for a larger range of eigenvalues in figure 6.6 where we plot 
log |^Af,exact " ^Af,WKB| as function of N. 
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Qualitative behaviour of the WKB expansion 
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Figure 6.5.: Behaviour of the WKB expansion of the few lowest eigenvalues E^. The blue 



m=o ^rn ^"^ i" equation (6.29). The green line 



points show the (asymptotic) sum J2 

shows the full expression t{6); it is independent of M to the expected numerical 



accuracy of the correction term (6.30). The red line shows the exact eigenvalue, 
evaluated numerically by our very-high-precision routine. The Borel corrected 
WKB series does not converge towards the exact eigenvalue. 

We find empirically that 

-E^Ar,exact > -E^Ar.WKB for all even N = 2M, 

with a difference which varies like e~'^^ with for A^ < 42, and approximately like e^'^{'^'^+^/'^) 
for iV > 42. Further 

^7V,exact < ^7V,WKB for odd iV = 2M + 1 < 9, 

also with a difference which varies like e~'^ with A^, and 

^iv,exact < ^JV,WKB for odd A^ = 4M + 3 > 11, 
Siv,exact > En^wkb for odd A^ = 4M + 1 > 13. 
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In these cases the difference behaves approximately like e' 



-7r(8+Ar/4) 



Table 6.1.: Behaviour of WKB expansion for Eq 



M 


-^0,WKB 


p(M-l) p(M) 
-^0,WKB + -^O,corr 


^0, exact 


1 
2 
3 

4 


0.867145 326484 821 
0.989 821295 452 906 
0.940 878 506 803 713 
0.842 885181871221 


0.949 048 242147079 
0.949 048 242 213 528 
0.949 048 245 949142 
0.949 048 880 595 005 


1.060 362 090484183 
1.060 362 090484183 
1.060 362 090484183 
1.060 362 090484183 



Behaviour of the WKB series for the lowest eigenvalue Eq. The first column shows the results 
of summing the first M terms of the WKB-series. The second column the result after the 



correction term (6.30) from Borel resummation has been added. The last column shows the 
eigenvalue computed numerically to very high precision. 





Table 6.2.: Be 


haviour of WKB expansion for Ei 


M 


-^l.WKB 


p(M-l) p(M) 
-^1,WKB ^ -^l,corr 


-C/l, exact 


1 


3.751919 923 550433 


3.808 235 541533 203 


3.799 673 029 801394 


2 


3.810 896 378 060 855 


3.808 235 541533 340 


3.799 673 029 801394 


3 


3.808 282 018 212 746 


3.808 235 541531506 


3.799 673 029 801394 


4 


3.807 700409 855 639 


3.808 235 541531468 


3.799 673 029 801394 


5 


3.808 311380 649 850 


3.808 235 541532 831 


3.799 673 029 801394 


6 


3.808 972 737814 702 


3.808 235 541513 511 


3.799 673 029 801394 


7 


3.807487485 686 370 


3.808 235 542141864 


3.799 673 029 801394 


8 


3.803436 692 708 719 


3.808 235 536164 707 


3.799 673 029 801394 



Behaviour of the WKB series for the eigenvalue Ei. The first column shows the results of sum- 
ming the first M terms of the WKB-series. The second column the result after the correction 



term (6.30) from Borel resummation has been added. The last column shows the eigenvalue 
computed numerically to very high precision. 



It should be clear that the Dunham quantization formula (4.21 ) does not provide exact 



eigenvalues in this case. There are (leading) order correction terms which look intriguingly 
simple. They are manifestations of the fact that the WKB approximation is inexact, even 
when summed to arbitrarily high order. Starting with a WKB-solution which behaves like 



Q(z)-V4exp(^-i£dt^^Q(^+...) 
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no higher-order correction will provide a contribution which changes the sign of the square 
root, i.e., provides a solution which behave like 



Q(z)-i/^exp0£dt/^Q(^ + •••). 



However, both behaviours are usually present in the exact solution. In asymptotic analysis 
they are said to emerge when Stokes lines are crossed. 



Logarithmic difference between exact and WKB eigenvalues En 
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Figure 6.6.: Difference between the exact eigenvalues E]\j^exact (computed numerically to very 
high precision) and the WKB eigenvalues £^Ar,wKB. computed using either the 
optimal asymptotic approximation (OAA) or adding the correction integral from 
Borel resummation (Borel). The results of these two methods cannot be distin- 



guished in the figure when A^ > 1. The later is found from equation (6.11), with 



t{5) computed from equation (6.29) for a range of Af -values around 2.2N. The 
result varies little with M, as indicated by the plotted standard deviation (t{En). 
Hence, the difference between -Ejv exact and Ef^\j\i\<^^ is much larger than the un- 



certainty in £^Ar,wKB due to numerical evaluation of the integral (6.30), although 



exponentially small as function of N. The correction terms look quite simple, 
with an interesting difference between the even and odd eigenvalues. 



Bibliography 

[1 

[2 



Alexander J. Yee, Shigeru Kondo et. al., http://www.numberworld.org/ 
niisc_runs/pi-10t/details .html (2011) 



The Royal Swedish Academy of Sciences, http://www.nobelprize.org/ 

nobel_prizes/physics/laureates/2005/ 

advanced-physicsprize2005.pdf 

[3] J. Beringer et. al. (Particle Data Group) The Review of Particle Physics, 
Phys. Rev. D 86, 010001 (2012) 

[4] T. Aoyama, M. Hayakawa, T. Kinoshita, and M. Nio, Tenth-Order QED Contribution 
to the Electron g — 2 and an Improved Value of the Fine Structure Constant, Phys. 
Rev. Lett. 109, 111807 (2012) 

[5] T. Aoyama, M. Hayakawa, T. Kinoshita, and M. Nio, Complete Tenth-Order QED 
contribution to the Muon g-2, Phys. Rev. Lett. 109, 111808 (2012) 

[6] L. Fuchs, Zur Theorie der linearen Differentialgleichungen mit verdnderlichen Coeffi- 
cienten. Journal fiir die reine und angewandte Mathematik 66, 121 (1866) 

[7] F.G. Frobenius, Uber die Integration der linearen Differentialgleichungen durch Rei- 
hen. Journal fiir die reine und angewandte Mathematik, 76, 214 (1873) 

[8] Philip M. Morse and Herman Feshbach, Methods of Theoretical Physics, Part I, 508- 
515, McGraw-Hill (1953) 

[9] Philip M. Morse and Herman Feshbach, Methods of Theoretical Physics, Part I, 663, 
McGraw-Hill (1953) 

[10] A.D. Fokker, Die mittlere Energie rotierender elektrischer Dipole im Strahlungsfeld, 
348, Ann. Phys. 810-820 (1914) 

[11] M. Planck, Sitz.ber (1917) 

[12] B. Haible and R.B. Kreckel, CLN - Class Library for Numbers, 
http : //www . ginac . de/CLN/ 

[13] T. Granlund and collaborators, CMP - The CNU Multiple Precision Arithmetic Li- 
brary, http://gmplib.org/ 



57 



58 Bibliography 

[14 
[15 



A. Schonhage and V. Strassen, Schnelle Multiplikasjon grofier Zahlen, Computing 7, 
281-292 (1971) 



J. Zinn-Justin and U.D. Jentschura, Multi-Instantons and Exact Results I: Conjec- 
tures, WKB Expansions, and Instanton Interactions, Annals of Physics 313, 197-267 
(2004) arXiv:quant-ph/0501136 

[16] A. Mushtaq, A. Kv^rno and K. Olaussen, Systematic Improvements of Splitting Meth- 
ods for the Hamilton Equations, Proceedings of The World Congress on Engineering 
2012 Vol I, WCE 2012, July 4-6, 2012, London, U.K., 247-251. arXiv: 1204.4117 

[17] Niels Bohr, On the Constitution of Atoms and Molecules, Philosophical Magazine 
Series 6, 26, 1-25 (1913) 

[18] Niels Bohr, On the Constitution of Atoms and Molecules Part II, Philosophical Mag- 
azine Series 6, 26, 476-502 (1913) 

[19] William Wilson, The quantum theory of radiation and line spectra. Philosophical Mag- 
azine Series 6, 29, 795-802 (1915) 

[20] Arnold Sommerfeld, Zur Quantentheorie der Spectrallinien, Annalen der Physik 51, 
1 (1916) 

[21] H. Jeffreys, On certain approximate solutions of linear differential equations of the 
second order. Proceedings of the London Mathematical Society 23, 428-436 (1924) 

[22] G. Wenzel, Eine Verallgemeinerung der Quantenhedingungen fAijr die Zwecke der 
Wellenmechanik, Zeitschrift fiir Physik 38, 518-529 (1926) 

[23] H.A. Kramers, Wellenmechanik und halbzdhlige Quantisierung, Zeitschrift fiir Physik 
39, 828-840 (1926) 

[24] L. Brillouin, La mecanique ondulatoire de Schrodinger: une methode generale de res- 
olution par approximations successives, Comptes Rendus de I'Academie des Sciences 
183, 24-26 (1926) 

[25] L.I. Schiff, Quantum Mechanics, Third Edition, section 34, McGraw-Hill (1968) 

[26] H. Kroemer, Quantum Mechanics: for engineering, materials science, and applied 
physics. Chapter 6, Prentice Hall (1994) 

[27] CM Bender and S.A. Orszag, Advanced Mathematical Methods for Scientists and 
Engineers, Chapter 10, McGraw-Hill (1978) 

[28] J.L. Dunham, The Wentzel- Brillouin- Kramers Method of Solving the Wave Equation, 
Phys. Rev. 41, 713 (1932) 

[29] CM. Bender, K. Olaussen, and P.S. Wang, Numerological analysis of the WKB ap- 
proximation in large order. Physical Review D16, 1740-1748 (1977) 



59 



[30] Paul B. Bailey, Exact Quantizatization Rules for the OneDimensional Schrodinger 
Equation with Turning Points, Journal of Mathematical Physics 5, 1293-1297 (1964) 

[31] C. Rosenzweig and J.B. Krieger, Exact Quantization Conditions, Journal of Mathe- 
matical Physics 9, 849-860 (1968) 

[32] Philip M. Morse, Diatomic molecules according to the wave mechanics. II. Vibrational 
levels, Phys. Rev. 34, 57-64 (1929) 

[33] G. Poschl and E. Teller, Bemerkungen zur Quantenmechanik des anharmonischen 
Ozillators, Zeitschrift fiir Physik 83, 143-151 (1933) 

[34] Serge Lang, Linear Algebra Third Edition, 262-264, Springer- Verlag (1987) 

[35] Asif Mushtaq, Amna Noreen, Kare Olaussen, Ingjald 0verb0, Very-high-precision 
solutions of a class of Schrodinger type equations, Computer Physics Communications 
182, 1810-1813 (2011); arXiv: 1008.0834 

[36] Amna Noreen and Kare Olaussen, Very-high-precision normalized eigenfunctions for 
a class of Schrodinger type equations. Proceedings of World Academy of Science, En- 
gineering and Technology 76, 831-836 (2011); arXiv: 1105. 1460 

[37] Amna Noreen and Kare Olaussen, High precision series solution of differential equa- 
tions: Ordinary and regular singular point of second order ODEs, Computer Physics 
Communications 183, 2291-2297 (2011); arXiv: 1205.2226 

[38] Amna Noreen and Kare Olaussen, Estimating Coefficients of Frobenius Series by Leg- 
endre Transform and WKB Approximation Proceedings of the World Congress of 
Engineering 2012 Vol II WCE 2012, London U.K. 789-791; arXiv: 1205.2221 

[39] Amna Noreen and Kare Olaussen, Cenerating Very-High-Precision Frobenius Series 
with Apriori Estimates of Coefficients, lAENG International Journal of Computer 
Science 39, 386-393 (2012); arXiv: 1209.6237 

[40] Amna Noreen and Kare Olaussen, Quantum loop expansion to high orders, extended 
Borel summation, and comparison with exact results, arXiv: 1209.6242 



